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ABSTRACT 

Recently  Kundu( 1993a)  has  proposed  a  non-linear 
eigenvalue  method  for  finding  the  maximum  likelihood 
estimators  (MLE)  of  the  parameters  of  undamped 
exponential  signals.  It  is  known  to  perform  better 
than  the  previously  existing  methods,  like  FBLP  of  Tufts 
and  Kumaresan  (1982)  or  Pisarenko's  method  (Pisarenko; 
1972),  in  the  sense  of  lower  mean  squared  errors.  The 
solution  in  general  depends  on  the  roots  of  a 

polynomial  equation.  It  is  observed  that  the 

coefficients  of  the  polynomial  exhibit  a  certain 
symmetry.  Since  it  is  known  (Crowder;  1984)  that  the 
MLE  with  constraints  is  more  efficient  than  the 

unconstrained  MLE,  modified  maximum  likelihood  method 
has  been  suggested  to  estimate  the  parameters  under 
these  symmetric  constraints.  It  is  observed  in  the 
simulation  study  that  the  mean  squared  errors  of  the 
constrained  MLE  are  closer  to  the  Cramer-Rao  lower 
bound  than  the  ordinary  MLE  in  almost  all  the 
situations. 

Key  Words  and  Phrases:  Constrained  MLE, 

Sinusoids,  Prony's  Algorithm,  EVLP,  FBLP. 


Multiple 


1. 


INTRODUCTION 


In  signal  processing  detecting  the  number  of 
signals  and  estimating  the  parameters  of  the  signals 
are  very  important  problems.  In  this  article,  we 
consider  the  estimation  of  the  following  exponentially 
undamped  sinusoidal  signal  with  additive  noise,  i.e. 


M 

Yn  =  “k  ®  ^ . “ 

Here  are  unknown  complex  numbers,  e  [0,2n)  are 

unknown  radian  frequencies  and  i  ='  The  are 

independent  identically  distributed  (i.i.d.)  complex 
valued  error  random  variables  such  that  E(€^)  =  0  and 

E  j  I  =  O'  .  oi  ' ' '  ’  r  ^  ~  ^^1  ^  *  *  *  ure 

unknown  and  the  are  assumed  to  be  distinct.  The 

problem  is  to  estimate  the  unknown  parameters,  given  a 
sample  of  size  N.  Here  we  assume  that  M,  the  number  of 
signals,  is  known. 

Furthermore  consider  the  following  non-linear  time 
series  regression  model; 


y 


n 


L 

J]  [Aj^  cos(coj^n)  +  sin(tOj^n)  ]  + 


n  =  1, . . . ,N  (1.2) 


Here  and  Bj^are  real  amplitudes  and  are 
frequencies  and  are  i.i.d.  real  valued  random 
variables  with  mean  zero  and  finite  variance  . 
Observe  that  the  model  (1.2)  can  be  written  as  a 
special  case  of  the  model  (1.1). 

The  problem  of  estimating  the  parameters  of 
complex  sinusoids  in  noise  has  received  considerable 
attention  in  the  past  several  years^  /see  for  example 
Kay  and  Marple  (;L981)and  Kay;  (l988jchapter  13^  in  the 


2 


/  J 


Electrical  Engineering  literature  where  as  a  reasonable 
amount  of  work  has  been  done  starting  with  the  work  of 
Hannan  (1971)  and  Walker  (1971), in  the  Statistics 
literature  also^  (see  for  example  Rice  and  Rosenblatt; 
^9 8^ and  Kundu  1993b ^f or  estimating  real  sinusoids. 
Under  the  normality  assumptions  on  the  error  term,  the 
maximum  likelihood  estimators  (MLE)  of  u  and  a  are  same 

as  the  non-linear  least  squares  estimators  obtained  by 
minimizing 


N 

RnCS/®)  =  L 
~  ~  n=l 


n 


M 

-  E 

k=l 


in©. 


^k 


(1.3) 


with  respect  to  ©  =  (6  ^,...,9^)  and  a  =  (a^ , . . .  ,3^^) . 

The  solution  obtained  from  (1-3)  is  a  consistent 
estimator  of  co  and  a  (Mitra  and  Kundu;  1993  or  Rao  and 

Zhao  (1993)).  Unfortunately  it  is  known  that  (Smyth; 
1985,  Kundu;  1989)  the  general  purpose  algorithms  such 
as  Gauss  Newton,  Newton  Rapson  or  their  variants  often 
have  great  difficulty  in  converging  to  the  optimum 
solution.  Recently  Kundu  (1993a)  gives  an  efficient 
iterative  procedure  to  obtain  the  MLE  of  a>,  the 

asymptotic  stability  of  the  procedure  can  be  found  in 
Kundu  (1994). 

In  estimating  the  frequencies,  the  methods  like 
the  Pisarenko's  method  (Pisarenko;  1972),  Forward 
Backward  Linear  Prediction  (FBLP)  of  Tufts  and 
Kumaresan  (1982),  EquiVariance  Linear  Prediction  (EVLP) 
discussed  in  Bai,  Krishnaiah  and  Zhao  (1986)  or  -Rao 
(1988),  are  based  on  the  classical  methods  of  Prony 
(1795).  The  basic  MLE  Prony  method  goes  back  to  Osborne 
(1975).  That  method  was  further  developed  and 

investigated  by  Osborne  and  Smyth  (1991,  1994)  and  Kahn 
et.  al.  (1992).  Bai,  Rao  and  Chow  (1989)  refine  the 
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EVLP  estimators  to  produce  efficient  estimators  of  the 
frequencies.  The  solution  in  general  depends  on  the 
roots  of  a  polynomial  equation.  It  is  observed  that 
the  coefficients  of  the  polynomial  equation  satisfy 
some  symmetric  constraints.  Since  it  is  known 
(Crowder;  1984)  that  the  constrained  MLE  has  lower  mean 
squared  error  then  the  unconstrained  MLE,  we  introduce 
constrained  MLE  for  this  problem  and  compare  it  with 
the  existing  methods  using  simulation.  Observe  that 
once  we  obtain  the  non-linear  parameters  the  linear 
parameters  can  be  obtained  using  the  simple  linear 
regression  technique  |^Kundu(^  1994),. 

The  rest  of  the  paper  is  organized  as  follows,  we 
provide  the  Pisarenko's  method  in  Chapter  2.  The 
Modified  Prony  algorithm  is  given  in  Chapter  3,  Chapter 
4  contains  the  Constrained  maximum  likelihood 
procedure.  Some  Numerical  results  are  provided  in 
Chapter  5  and  we  draw  conclusions  from  our  results  in 
Chapter  6.  In  the  Appendix  we  provide  some  theoretical 
justification  of  the  convergence  of  the  numerical 
procedures  under  the  assumptions  of  independent  and 
identically  distributed  error  random  variables. 

2.  PISARENKO'S  METHOD 

Prony  (1795)  suggested  a  method  of  solving  the 
non-linear  estimation  problem.  Many  standard  texts  on 
numerical  methods  outline  this  algorithm  (Barrodale 
and  Olesky  (1981),  Froberg  (1969)).  The  algorithm  can 
be  extended  to  the  noisy  case  as  follows  ; 

Suppose  that  the  vector  g  =  (g^, . . . .  ,gj^)'  is 

such  that 

M  ^ 

go  +  z  +  ...  +  g„  z  =  g„^n(z  -  e  )  (2.1) 
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then  for  any  n  >  M+l 


k=0®'"  k=0®'' 


(2.2) 


where  the  right  hand  side  of  (2.2)  is  a  function  of 
error  only.  The  coefficients  gj  are  estimated  by 
minimizing 


N 


M 


E  I  E  gv  ^n-kl 

n=M+l  k=0  ^  ^  ^ 


(2.3) 


subject  to  the  condition  |g|  =  1.  Such  a  method  of 

estimation  is  known  as  the  Pisarenko's  method^  (also 
named  as  EVLP  method  by  Bai,  Krishnaiah  and  Zhao/ 
^986) . 

Now  write 


T  = 


'N-M 


'M+1 


■N 


(2.4) 


and  R  =  T  T,  here  '*  '  denotes  the  complex  conjugate 
transpose  of  a  matrix  or  of  a  vector. 


Therefore  the 


estimator  of  g  can  be  obtained  by  minimizing  g  Rg  such 

*  ~  A  _  ~  ~ 

that  g  g  =  1.  If  g  is  the  estimator  of  g  so  obtained 


construct  the  polynomial  equation 

A  A  ^  M 

Q  +  a,  z  + . . .  +  a„  z  =0 
^o  ^1 


(2.5) 


obtain  solutions  in  the  form  of 


A 

-i<o. 


-ICi). 


Pi  e 


./  Pjt  ® 


'M 


with  Pi  >  0 


for  i  =  1,2, .. .  M. 


(2.6) 


and  take  aji,...,a>j^  as  estimators  of  , . . .  It  is 
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shown  in  Bai,  Krishnaiah  and  Zhao  (1986)  that  a>  is  a 

consistent  estimate  of  co  with  a  convergence  rate  of 

0  It  is  observed  that  t^ao  and  Zhao;  (1993^ 

P  A 

Mitra  and  KundUj-^^  1993^  co  obtained  by  minimizing  (1.2) 

is  a  consistent  estimator  of  a>  with  a  convergence  rate 

of  Op(N~^^^).  On  the  other  hand  the  FBLP  of  Tufts  and 
Kumaresan  (1982)  may  not  provide  consistent  estimator 
of  CO  (Rao;  1988) . 

3 .  MODIFIED  PRONY  ALGORITHM 


A  modification  to  Prony's  algorithm  have  been 
considered  by  Marple  (1979,1987),  using  the  fact  that 
the  roots  of  the  polynomial  equation  (2.5)  should  be  of 
unit  modulus.  Consider  the  polynomial  equation 

p(z)  =  go  +  g^z  +  •••  +  gj^  =  0  (3.1) 

with  Zj^,  k  =  1,....M  as  its  roots.  Since  the  Zj^'s 
are  of  unit  modulus,  we  have  z~^  =  Therefore  the 
polynomial  P(z)  and 


R(2)  =  g^i  +  %-!  z  +...+  g^  z 


M 


(3.2) 


have  the  same  roots.  Hence 

9k  % 


^M- 


M-K  9o 


or 


’K 


'O 


^M 


%-k!  - 
^o 


- 1 
J 


1/2 


K  =  0, . . . ,M 


K  =  0,  .  .  .M 


(3.3) 


If  we  define 
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‘’k  "  % 


,1/2 


1£ 

% 


K  =  0, . . .M 


(3.4) 


then  it  is  easy  to  see  that 

Thus  (2-2)  can  be  written  as;  if  Yj^  satisfies  (1.1), 
then  there  exists  a  vector  b  =  (t>Q/  •  •  •/  bj^)  such 

that  bj^  =  ^M-k  ^  =  0 , .  . .  ,M  and 


M 

2  b 
k=0 


k 


^n-k 


M 

Z  ] 
k=0 


k  n-k 


(3.5) 


for  any  n  >  M+1. 


4 .  CONSTRAINED  MLE 

Bresler  and  Macovski  (1986)  and  Kumaresan,  Scharf 
and  Shaw  (1986)  considered  the  maximum  likelihood 
approach  using  similar  conditions  on^-the  coefficients. 
They  transformed  the  problem  to  a  constrained  quadratic 
minimization  problem  and  solved  it  iteratively  using 
some  standard  computer  packages.  No  proof  of 

convergence  was  provided  in  these  papers.  In  this 
section  we  modify  the  MLE  procedure  of  Kundu  (1993a) 
using  similar  conditions  on  the  coefficients.  We 
transform  the  problem  to  a  non-linear  eigenvalue 
problem  and  solve  it  iteratively. 

The  model  (1.1)  can  be  written  in  the 
following  matrix  form; 

Y  =  A(a5)oi  +  E  (4.1) 

where  y'^  =  (Y^^  /  •  •  • /Yj^)  /  =  (a^^ ,  .  .  .  ,0^^) ,  = 

~  T  ~  ~ 

(<o^ ,  .  .  -  ,Ci>j^)  and  E  =  (e^,...,e^).  A(<»))  is  a  NxM 

^“gP 

matrix  and  its  (p,q)th  element  ise  ^,p=l,...,N,  q 
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=  1, . . . ,M. 

Therefore  (1.3)  can  be  written  as 

Rj^(a,  oj)  =  |Y  -  A(aj)oi|^  (4.2) 

Here  the  linear  parameters  ot  are  separable  from  the 
non-linear  parameters  co.  For  a  fixed  value  of  cj,  the 
minimization  of  with  respect  to  oi  is  a  simple  linear 

A 

regression  problem.  The  solution  of  a.(u>)  is 

a(a>)  =  (A*(a>)  A(a>))‘^  A*(co)Y  (4.3) 

Substituting  back  a  (to)  in  (4.2)  and  denoting  Q(w)  = 
Rjj(a,a))  we  obtain 

Q(to)  =  Y*(I  -  P^)Y  (4,4) 

where  =  A(A*A)~^A*  is  the  projection  operator  in 

the  range  space  spanned  by  the  columns  of  the  matrix  A. 
Therefore  the  MLE  of  <o  can  be  obtained  by  minimizing 

(4.4)  ,  w.r.t.  CO. 

Consider  the  NxN-M  matrix  X  of  the  following  form: 


(4.5) 


X  = 


0 

0 


o 


0 

0 

0 

* 


From  (2.2)  it  follows  X  =  0,  where 
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(Ey^, .  . .  ,Eyj^)  .  This  implies  X  ACoj)  =  0,  i.e.  the 

columns  of  A  are  orthogonal  to  the  rows  of  X  and  I 

=  X(X*X)~^X*.  Therefore  minimizing  Y*(I-P^)Y 

.  .  .  .  ~  *  ~ 

with  respect  to  co  is  equivalent  to  minimizing  Y  P^^  Y 
with  respect  to  g  =  (g^/ .  •  •  ,gj^)  • 

A 

Let  g  be  the  vector  that  minimizes 

Y*X(X*X)”^X*Y.  Then  the  MLE  of  to  can  be  obtained 

by  solving  the  polynomial  equation  as  described  in 
(2.5)  and  (2.6)  ^see  Kundu;  (l993a) .  Similarly  the 

constrained  MLE  of  co  can  be  obtained  by  minimizing 


"k  k  —1  k 

Q(g)  =  Y  X(X  X)  X  Y 


(4.6) 


such  that  g*g  =  1  and  g^^  =  ;k  =  0,...,M. 

Then  the  estimator  of  w  can  be  obtained  by  solving  the 


polynomial  equation  (2.5).  Observe  that  the  roots  of 
the  polynomial  equation  in  the  constrained  case  will  be 
of  the  form 


A 


e 


(4.7) 


instead  of  (2.6). 

The  constrained  minimization  can  be  done  as 
follows.  The  matrix  X  in  the  constrained  case  can  be 
written  as  X  =  [x^ , . . . ,  where  x^^  is  of  the  form  [ 

0,g,0  ].  Let  U  and  V  denote  the  real  and  imaginary 

parts  of  X  respectively.  Assuming  M  is  odd  for 
brevity  and  k  =  0,  1,...,  (M-l)/2. 
Then  X  can  be  written  as 


X 


(M-l)/2 

=  E 

a=0 


(c  U 
a  a 


+  id  V  ) 
a  a 


(4.8) 
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where  U  ,  V  ;  a  =  0,1,...,  are 

a.  a.  '  '  '  2 


NXN-M 


matrices 


with  entriesO  and  1.  The  minimisation  of  Q(g)  = 

^  %-k' 


Y*X(X*X)”^X*Y  such  that  g*g  =  1  and  g,  =  g^-v'  ^  = 


0,...,M  is  equivalent  to  minimize 

Q(c,  d)  +  \  (c'^c  +  d^d  -  1) 


(4.9) 


with  respect  to  c,  d  and  X.  Here 


we 


denote  c 


7  _ 


/  M-l.  ,T 

-7-^'  ^  " 


d. 


(d^,...,  Q{c,  d) 


Y*X(X*X)  ^X*Y,  where  X  is  of  the  form  (4.8)  and  X  is 

the  Lagrange  multiplier.  If  we  differentiate  (4.9) 
with  respect  to  c  and  d,  we  obtain  a  matrix  equation  of 

the  form 


D(c^,d^) 


2X 


=  0 


(4.10) 


Here  D  is  a  M+lxM+1  matrix.  The  matrix  D  can  be  written 
in  a  partitioned  form  as  : 


D 


B 


B 

C 


(4.11) 


M+i  M+1 

where  A,B,C  are  all  — ^  x  — ^  matrices.  We  have 


=  Y*Uj(X*X)  ^U*Y  -  Y*X(X*X)  ^(UjUj^+  U^U  ^ )  ( X*X)  "  Vy 
+  Y*Uj^(X*X)“^UjY  j,k  =  0,1,...  (4.12) 

Bjk  =  -iY*Uj(X*X)"^V*Y-iY*X(X*X)"^(U^Vj^-V^Uj)(X*X)"VY 

.  *  *  -1  T  .  M-1 

+  lY  Vj^(X  X)  UjY  3,k  =  0,1,...  (4.13) 
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C.y.  =  Y*Vj(X*X)  ^V*y-Y*X(X*X)  +  V^Vj)(X*X)  ^X*Y 


*  *  -1  T 

+  Y  V,  (X  X)  vtY 
■K  J 


j,k  —  0,1,... 


M-1 


(4.14) 


T  T 

Premultiplying  the  left  hand  side  of  (4.10)  by(c  ,  d  ), 


we  obtain 


[c'^,  d*^]  D(c^,d^) 


+  2X  =0 


(4.15) 


It  can  be  easily  seen  that 


Io\  d'^l  0(0’^,^^) 


=  0 


(4.16) 


which  implies  X  =  0.  Therefore  the  constrained 

minimization  of  Q(g)  is  equivalent  to  solving  a  matrix 


equation  of  the  form 
c 


D(c^,d^) 


=  0 


T  ^  T 

s.t.  cc+dd=l 


(4.17) 


T 

similarly  when  M  is  even,  c  =  (c^, 


and  D  can  be  partitioned  as 


D  = 


B 

C  . 


/ 


(4.18) 


M  ,  M  ,  ,  M  .  ,  MM  M 

Here  A ,  B ,  C  are  1x^+1,  ^  "il  x  2 ;  2^2^ 

matrices  respectively.  The  (j,k)th  element  of  A,  B, 

or  C  is  same  as  before  with  appropriate  ranges.  Note 

that  D  is  a  real  symmetric  matrix  in  both  the  cases, 

and  we  need  to  solve  a  matrix  equation  of  the  form 


D(x)x  =  0  s.t.  Ilxll  =  1  (4.19) 
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This  is  a  nonlinear  eigenvalue  problem.  x 

satisfying  (4.19)  should  be  an  eigenvector 
corresponding  to  the  zero  eigenvalue  of  the  matrix 
D(x) .  We  suggest  the  following  iterative  technique 


similar  to  that  of  Osborne  (1975)  and  Kundu  (1993a)  to 
solve  (4.19)  : 


(D(x^)  -  \^‘''^I)x^'''^  =  0 


(4.20) 


where  is  the  eigenvalue  closest  to  zero  of  D(x  ) 

and  x^^^  is  the  corresponding  normalized  eigenvector. 

~  J  1«  k+1 

The  iterative  process  should  be  stopped  when  x  is 

small  corresponding  to  11  Dll  . 

The  algorithm  has  the  following  form  : 

Step  1:  Set  an  initial  value  x^  and  normalize 
it,  i.e. 

;  i  =  1. 


Step  2:  Calculate  matrix  D(x^) 

Step  3:  Find  the  eigenvalue  of  D(x^) 

closest  to  zero  and  normalize  the 
corresponding  eigenvector  x^^^ 

Step  4:  Test  the  convergence  by  checking  if 

|X^'‘'^|  <  e  II  Dll 

Step  5:  If  the  test  in  4  fails,  i:=  i+1  and  go  to 
Step  2 . 
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5 .  NUMERICAL  EXPERIMENTS : 


We  performed  some  numerical  experiments  mainly  to 
compare  the  performances  of  the  different  methods  for 
finite  samples.  All  these  simulations  are  done  using 
the  IMSL  (1984)  random  deviate  generator.  We  consider 
the  following  two  models: 


Model  1. 


^  „  l.Oin  ,  o  n  2.0iin  , 
y  =  2 . 0  e  +  3 . 0  e  +  e. 

n 


n 


n  =  1 , . . . ,N. 

Model  2 . 

y^  =  1.0  cos (2. On)  +  1.0  sin (2. On)  + 


(5.1) 


(5.2) 


For  Model  1,  e  are  i.i.d.  complex  valued  normal  random 

variables  with  mean  zero  and  standard  deviation  a  for 

both  the  real  and  imaginary  parts.  The  real  and 

imaginary  parts  are  taken  to  be  independently 

distributed.  For  Model  2  e  are  i.i.d.  real  valued 

n 

random  variables  with  mean  zero  and  standard  deviation 
a.  For  different  values  of  N  and  a,  one  thousand 
different  data  sets  were  generated.  Numerical  results 
are  observed  for  sample  sizes  N  =  20,30,40,50  and  a  = 
.01,.!  and  .5.  Observe  that  our  method  can  be  applied 
even  if  the  errors  are  stationary  random  variables. 
For  Model  2,  we  also  consider  when  are  moving 

average  error  of  the  form; 


where  e._''s  are  i.i.d.  normal  random  variables  with  mean 


n 


40 


zero  and  standard  deviation  o'  =  .1  and  N  =  20,  30, 
and  50 . 

For  each  data  set,  we  computed  the  ordinary  MLE 
(MLE)  as  described  in  Kundu  (1993a),  constrained  MLE  as 
in  Section  4  (CMLE)  and  as  in  Bresler  and  Macovski 
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(CMLEBM),  the  FBLP  as  in  Tufts  and  Kumaresan  (1982)  and 
the  modified  EVLP  (EVLP)  as  in  Bai,  Rao  and  Chow 
(1989).  For  the  MLE  and  CMLE  we  use  s  =  10~**,  the 
ordinary  EVLP  estimator  as  the  starting  value  and  ||D| 
to  be  the  largest  eigenvalue  of  the  matrix  D-  For 
CMLEBM,  we  use  the  same  initial  value  as  MLE  or  CMLE 
and  7  steps  iterations  in  all  the  cases.  Observe  that 
the  stopping  rule  of  MLE  or  CMLE  is  quite  different 
than  that  of  CMLEBM.  It  is  observed  that  for  MLE  or 
CMLE  the  iteration  converged  in  7  steps  in  all  the 
cases  for  both  the  models,  with  the  above  stopping 
criterion.  Therefore  to  have  a  reasonable  comparison 
of  MLE  or  CMLE  with  CMLEBM  we  use  7  steps  iterations 
in  all  the  cases  in  both  the  models  for  CMLEBM  also. 
The  mean  estimates  of  the  frequencies  and  their  average 
mean  squared  errors  (MSE)  over  one  thousand 
replications  are  reported. 

Observe  that  we  write  (5.2)  in  the  form  of  (1.1) 
with  2.0  and  =  -2.0.  We  estimated  and  by 
the  different  methods  as  described  above.  In  all  our 
simulations  for  all  the  methods  the  two  roots  of  (2.5) 
came  in  the  complex  conjugate  form,  we  obtained  the 
estimator  of  oj  from  that.  We  also  obtain  the 
estimators  of  the  linear  parameters  from  (4.3)  in  all 
the  cases. 

We  reported  all  the  results  in  Table  1-4  for  Model 
1  and  in  Table  5  for  Model  2  when  the  errors  are  i.i.d 
and  when  the  errors  are  from  (5.3),  the  results  are 
reported  in  Table  6 .  In  each  table  we  also  put  the 
Cramer-Rao  lower  bounds  (CRLB)  for  comparison. 
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Table  la 
Model  1 


N  = 

20,  a  =  . 

01 

Method 

\  ^ 

(0 

a 

a 

'  1 

z 

1 

z 

PARAMETERS 

1.00000 

2.00000 

2.00000 

3.00000 

(CRLB) 

( .194E-3) 

( .129E-3) 

( .224E-2) 

( .224E-2) 

MLE 

1.00007 

1.99989 

2.00025 

3.00010 

( .275E-3) 

( .201E-3) 

( .246E-2) 

( .247E-2) 

CMLE 

1.00000 

1.99997 

2.00012 

3.00005 

(.215E-3) 

( .146E-3) 

( .235E-2) 

( .242E-2) 

CMLEBM 

1.00001 

1.99997 

2.00010 

3.00005 

( .223E-3) 

(  .152E-3) 

( .240E-2) 

( .242E-2) 

FBLP 

1.00014 

2.00004 

3.99996 

5.50002 

( .560E-3) 

( .220E-3) 

( .287E-2) 

( .277E-2) 

EVLP 

0.99982 

2.00004 

2.00011 

3.00040 

( .586E-3) 

(  .217E-3) 

( .304E-2) 

( .280E-2) 

Table  lb 

N  = 

CO 

o 

"<« 

Q 

II 

• 

1 

Method 

\  w 

(0 

a. 

a. 

z 

1 

z 

PARAMETERS 

1.00000 

2.00000 

2.00000 

3.00000 

(CRLB) 

( .194E-2) 

(  .129E-2) 

( .224E-1) 

( .224E-1) 

MLE 

1.00009 

1.99983 

2.00044 

3.00007 

( .280E-2) 

( .197E-2) 

( .240E-1) 

( .251E-1) 

CMLE 

1.00000 

1.99992 

2.00010 

3.00005 

( .217E-2) 

(  .149E-2) 

( .230E-1) 

( .244E-1) 

CMLEBM 

1.00001 

1.99991 

2.00010 

3.00013 

( .225E-2) 

( .149E-2) 

( .230E-1) 

( .247E-1) 

FBLP 

2.00001 

3.49990 

3.99996 

5.50002 

( .531E-2) 

( .230E-2) 

( .311E-1) 

( .288E-1) 

EVLP 

2.00007 

3.49996 

3.99995 

5.50001 

( .521E-2) 

( .225E-2) 

( .312E-1) 

( .282E-1) 
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Table  Ic 

N  = 

20^  O'  • 

5 

Method 

\ 

Cl) 

2 

ot 

1 

a 

2 

PARAMETERS 

1.00000 

2.00000 

2.00000 

3.00000 

(CRLB) 

( .968E-2) 

( .645E-2) 

( .112E-0) 

( .112E-0) 

MLE 

1.00007 

1.99989 

2.00025 

3.00010 

( .176E-1) 

( .991E-2) 

( .201E-0) 

( .182E-0) 

CMLE 

1.00029 

1.99906 

1.98828 

2.99803 

( .108E-1) 

( .726E-2) 

( .123E-0) 

( .119E-0) 

CMLEBM 

1.00054 

1.99891 

1.98650 

2.99690 

( .113E-1) 

( .754E-2) 

( .132E-0) 

( .118E-0) 

FBLP 

0.98911 

2.01315 

1.85643 

2.88781 

( .391E-1) 

( .215E-1) 

( .340E-0) 

( .261E-0) 

EVLP 

0.98619 

1.99715 

1.76575 

2.88092 

.5.-1 

( .398E-1) 

( .229E-1) 

(  .324E-0) 

( .254E-0) 

Table  2a 

Model  1 

N  = 

30,  O'  =  .1 

01 

Method 

\ 

w 

z 

a. 

1 

a. 

z 

PARAMETERS 

1.00000 

2.00000 

2.00000 

3.00000 

(CRLB) 

( .105E-3) 

( .703E-4) 

( .182E-2) 

( .182E-2) 

MLE 

1.00007 

1.99989 

2.00025 

3.00010 

-:vl 

( .201E-3) 

( .lllE-3) 

( .202E-2) 

( .219E-2) 

■^1 

CMLE 

1.00001 

1.99999 

2.00016 

2.99916 

.3 

(  .140E-3) 

( .712E-4) 

( .193E-2) 

( .202E-2) 

/-■ 

CMLEBM 

1.00001 

1.99998 

2.00014 

2.99995 

( .146E-3) 

( .739E-4) 

( .193E-2) 

( .201E-2) 

FBLP 

0.99983 

2.00014 

1.99996 

3.00190 

( .299E-3) 

( .250E-3) 

( .270E-2) 

( .251E-2) 

•:*. 

EVLP 

0.99984 

1.99978 

2.00011 

3.00040 

-  •••r/ 

( .316E-3) 

( .276E-3) 

( .281E-2) 

( .277E-2) 
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Table  2b 
=  30 ,  a  ■=  . 


Method 

\  “ 

CO 

a 

'  1 

2 

1 

PARAMETERS 

1.00000 

2.00000 

2.00000 

(CRLB) 

( .105E-2) 

( .703E-3) 

( .182E-1) 

MLE 

1.00021 

2.00095 

2.00011 

( .201E-2) 

( .991E-3) 

( .231E-1) 

CMLE 

0.99989 

1.99989 

2.00009 

( .139E-2) 

( .722E-3) 

( .187E-1) 

CMLEBM 

1.00014 

1.99982 

2.00009 

( .145E-2) 

(  .741E-3) 

(  .191E-1) 

FBLP 

1.00011 

2.00019 

1.99563 

( .293E-2) 

( .199E-2) 

( .281E-1) 

EVLP 

0.99983 

1.99776 

2.01350 

( .310E-2) 

( .195E-2) 

(  .291E-1) 

Table  2c 

N  = 

:  30  ,  Cf  =  . 

5 

Method 

\  w 

iji 

a 

1 

2 

1 

PARAMETERS 

1.00000 

2.00000 

2.00000 

(CRLB) 

( .527E-2) 

( .351E-2) 

( .913E-1) 

MLE 

0.99987 

2.00012 

2.00011 

( .880E-2) 

( .481E-2) 

( .991E-1) 

CMLE 

1.00091 

1.99909 

2.04309 

(  .739.E-2) 

( .412E-2) 

( .920E-1) 

CMLEBM 

1.00094 

1.99913 

2.04420 

( .743E-2) 

( .416E-2) 

(  .928E-1) 

FBLP 

0.98763 

1.99763 

2.01341 

( .268E-1) 

( .175E-1) 

( .299E-0) 

EVLP 

0.98866 

1.98847 

1.81846 

( .276E-1) 

( .187E-1) 

( .306E-0) 

a 

2 

3.00000 
( .182E-1) 

2.99981 
( .243E-1) 

2.99891 
( .233E-1) 

2.99959 
( .211E-1) 

2.99856 
( .279E-1) 

2.98708 
( .277E-1) 


a 

2 

3.00000 
( .913E-1) 

3.00121 
( .972E-1) 

2.99037 
( .921E-1) 

2.99817 
(  .923E-1) 

2.97865 
(  .239E-0) 

2.78936 
(  .249E-0) 
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Table  3a 
Model  1 

N  =  40,  O'  =  -01 


Method 

\ 

CO 

2 

a 

1 

a 

2 

PARAMETERS 

1.00000 

2.00000 

2.00000 

3.00000 

1.', 

(CRLB) 

( .685E-4) 

( .456E-4) 

( .158E-2) 

( .158E-2) 

MLE 

1.00003 

1.99948 

2.00012 

2.99919 

( .891E-4) 

( .591E-4) 

( .189E-1) 

( .211E-1) 

'  J 

CMLE 

1.00001 

1.99998 

2.00091 

2.99960 

( .793E-4) 

( .495E-4) 

( .176E-1) 

( .173E-1) 

CMLEBM 

1.00001 

1.99998 

2.00019 

2.99987 

: 

( .806E-4) 

( .505E-4) 

( .176E-1) 

( .179E-1) 

FBLP 

1.00017 

2.00011 

1.99956 

2.99281 

( .248E-3) 

( .199E-3) 

( .199E-2) 

( .210E-2) 

m 

EVLP 

0.99899 

1.99981 

1.99962 

2.99981 

( .256E-3) 

( .196E-3) 

( .201E-2) 

( .202E-2) 

Table  3b 

N  = 

o 

Q 

II 

• 

1 

Method 

\ 

CO 

2 

a. 

1 

a. 

2 

PARAMETERS 

1.00000 

2.00000 

2.00000 

3.00000 

(CRLB) 

( .685E-3) 

( .456E-3) 

( .158E-1) 

( .158E-1) 

-r-'i 

MLE 

0.99979 

2.00091 

1.99911 

2.99891 

I- 

( .911E-3) 

( .685E-3) 

( .210E-1) 

( .206E-1) 

CMLE 

1.00012 

1.99981 

2.00889 

2.99698 

-tjs'-'a 

Y  ' 

( .793E-3) 

(  .496E-3) 

(  .175E-1) 

( .163E-1) 

CMLEBM 

1.00014 

1.99987 

2.00886 

2.99789 

( .801E-3) 

( .499E-3) 

( .181E-1) 

( .176E-1) 

FBLP 

0.99923 

2.00178 

1.99203 

3.00879 

y.Y^ 

( .269E-2) 

( .888E-3) 

( .271E-1) 

(  .259E-1) 

'"‘■r 

EVLP 

1.00076 

1.99786 

2.00171 

3.00932 

:-.Y. 

( .275E-2) 

( .885E-3) 

( .275E-1) 

(  .261E-1) 
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Table  3c 
N  =  40,  o-  = 


Method 

\ 

PARAMETERS 

1.00000 

(CRLB) 

( .342E-2) 

ii 

MLE 

0.99764 

( .423E-2) 

3 

CMLE 

1.00083 

( .395E-2) 

CMLEBM 

1.00089 
( .396E-2) 

.  r' 

FBLP 

1.09805 

''■'i 

( .218E-1) 

EVLP 

1.00165 

.  .-'Ki 

( .211E-1) 

-;-.g 

N  = 

Method 

\ 

-  ’ 

PARAMETERS 

1.00000 

■ 

(CRLB) 

( .489E-4) 

•  *1 

MLE 

0.99998 
( .601E-4) 

'-fi 

CMLE 

1.00000 
( .570E-4) 

CMLEBM 

0.99999 

(  .575E-4) 

FBLP 

1.00011 

( .158E-3) 

EVLP 

0.99994 

( .159E-3) 

hi 

a. 

a 

2 

1 

2 

2.00000 

2.00000 

3.00000 

( .228E-2) 

( .790E-1) 

( .790E-1) 

2.00120 

1.99910 

3.00019 

( .298E-2) 

( .899E-1) 

( .911E-1) 

1.99920 

2.03879 

2.98082 

( .248E-2) 

( .842E-1) 

( .817E-1) 

2.00011 

2.03563 

2.99914 

( .255E-2) 

( .851E-1) 

( .815E-1) 

2.00128 

2.02451 

3.01395 

( .936E-2) 

( .311E-0) 

( .231E-0) 

1.99765 

1.91876 

2.88976 

( .989E-2) 

( .299E-0) 

( .241E-0) 

Table  4a 

Model  1 

50 ,  O'  =  .  1 

01 

CO 

a 

a 

2 

1 

2 

2.00000 

2.00000 

3 . 00000 

( .326E-4) 

( .141E-2) 

( .141E-2) 

1.99988 

2.00009 

2.99999 

( .380E-4) 

( .181E-2) 

(  .178E-2) 

1.99999 

1.99998 

2.99991 

( .320E-4) 

( .171E-2) 

(  .140E-2) 

1.99999 

2.00001 

2.99990 

( .326E-4) 

( .176E-2) 

(  .141E-2) 

1.99920 

1.99931 

3.00016 

( .139E-3) 

( .189E-2) 

(  .189E-2) 

2.00010 

2.00078 

3.00001 

( .136E-3) 

( .198E-2) 

(  .198E-2) 
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Table  4b 

N  = 

50,  a  =  . 

1 

Method 

\ 

(0 

2 

a 

1 

a 

2 

PARAMETERS 

1.00000 

2.00000 

2.00000 

3.00000 

(CRLB) 

( .489E-3) 

( .326E-3) 

( .141E-1) 

( .141E-1) 

MLE 

1.00012 

2.00088 

1.99962 

3.00019 

Pi 

( .710E-3) 

(  .389E-3) 

( .191E-1) 

( .199E-1) 

CMLE 

0.99992 

1.99993 

2.00769 

2.99905 

( .670E-3) 

(  .335E-3) 

(  .167E-1) 

( .159E-1) 

CMLEBM 

1.00056 

2.00091 

2.00239 

3.00317 

( .675E-3) 

(  .340E-3) 

( .170E-1) 

( .169E-1) 

FBLP 

0.99918 

2.00098 

2.00270 

3.00349 

'W, 

( .238E-2) 

( .790E-3) 

( .262E-1) 

( .244E-1) 

EVLP 

1.00072 

1.99911 

2.00252 

3.00013 

( .244E-2) 

( .787E-3) 

( .251E-1) 

( .249E-1) 

Table  4c 

N  = 

:  50,  a  —  . 

5 

Method 

\ 

(0 

2 

a 

i 

a 

2 

:A 

PARAMETERS 

1.00000 

2.00000 

2.00000 

3.00000 

(CRLB) 

( .245E-2) 

( .163E-2) 

(  .707E-1) 

( .707E-1) 

MLE 

1.00569 

2.00253 

1.99987 

2.99911 

.•5 

( .375E-2) 

( .218E-2) 

( .849E-1) 

( .811E-1) 

-M 

CMLE 

0.99971 

1.99969 

2.03347 

2.99546 

( .327E-2) 

( .198E-2) 

( .830E-1) 

( .710E-1) 

CMLEBM 

0.99976 

2.00021 

1.99281 

2.99786 

' 

( .335E-2) 

( .205E-2) 

( .833E-1) 

( .707E-1) 

FBLP 

1.00198 

1.99875 

2.00320 

3.00237 

( .190E-1) 

( .788E-2) 

( .303E-0) 

( .225E-0) 

EVLP 

1.00880 

2.00210 

1.91642 

2.95185 

( .198E-1) 

( .779E-2) 

( .288E-0) 

(  .239E-0) 
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Table  5a 
Model  2 

a  =  .01,  CO  =  2.0 


Sample 

Methods 

Size\  20 

30 

40 

50 

CRLB 

2.00000 

2.00000 

2.00000 

2.00000 

( .387E-3) 

( .211E-3) 

( .137E-3) 

( .980E-4) 

MLE 

1.99989 

1.99996 

1.99995 

1.99989 

( .401E-3) 

( .311E-3) 

( .201E-3) 

( .123E-3) 

CMLE 

1.99996 

1.99997 

1.99998 

1.99999 

( .390E-3) 

( .219E-3) 

( .142E-3) 

( .988E-4) 

CMLEBM 

1.99998 

1.99999 

1.99998 

1.99998 

( .390E-3) 

( .220E-3) 

( .144E-3) 

( .995E-4) 

FBLP 

1.99987 

1.99990 

1.99989 

1.99995 

( .498E-3) 

( .437E-3) 

( .289E-3) 

( .241E-3) 

EVLP 

1.99988 

1.99989 

1.99999 

1.99998 

( .500E-3) 

( .431E-3) 

( .299E-3) 

( .238E-3) 

cr 

=  .  1  /  CO  =  2 

.0 

Sample 

Methods 

Size\  20 

30 

40 

50 

CRLB 

2.00000 

2.00000 

2.00000 

2.00000 

( .387E-2) 

( .211E-2) 

( .137E-2) 

( .980E-3) 

MLE 

1.99983 

1.99989 

1.99991 

2.00001 

( .403E-2) 

( .314E-2) 

( .199E-2) 

( .119E-2) 

CMLE 

1.99993 

2.00061 

1.99991 

2.00003 

( .388E-2) 

( .220E-2) 

( .141E-2) 

( .990E-3) 

CMLEBM 

2.00017 

1.99992 

1.99999 

1.99992 

( .388E-2) 

(  .224E-2) 

( .141E-2) 

( .994E-3) 

FBLP 

2.00011 

2.00032 

1.99994 

1.99988 

( .501E-2) 

(  .441E-2) 

( .286E-2) 

( .244E-2) 

EVLP 

2.00023 

1.99965 

1.99989 

1.99988 

(  .497E-2) 

( .434E-2) 

( .297E-2) 

( .240E-3) 
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1 


O’ 

• 

CM 

II 

3 

LT) 

• 

!1 

0 

Sample 

Size\  20 

30 

40 

50 

Methods 

CRLB 

2.00000 

2.00000 

2.00000 

2.00000 

( .194E-1) 

( .105E-1) 

( .685E-2) 

( .490E-2) 

MLE 

2.00011 

2.00004 

1.99995 

1.99997 

( .225E-1) 

(  .218E-1) 

( .789E-2) 

( .588E-2) 

CMLE 

2.00071 

1.99918 

1.99954 

2.00013 

■:  yj 

( .196E-1) 

( .109E-1) 

( .699E-2) 

( .498E-2) 

1 

i 

CMLEBM 

1.99901 

2.00018 

2.00001 

2.00071 

( .198E-1) 

(  .llOE-1) 

( .705E-2) 

( .500E-2) 

FBLP 

2.00088 

2.00098 

1.99934 

2.00108 

:i 

( .401E-1) 

(  .399E-1) 

( .303E-1) 

( .255E-1) 

EVLP 

1.99975 

1.99955 

2.00152 

1.99917 

( .402E-1) 

( .388E-1) 

( .298E-1) 

( .266E-1) 

Table  5b 

■  ! 

Model  2 

O' 

=  .01,  A  =  1 

.0 

Sample 

Size\  20 

30 

40 

50 

Methods 

‘ -  i 

CRLB 

1.00000 

1.00000 

1.00000 

1.00000 

( .500E-2) 

( .408E-2) 

( .354E-2) 

( .316E-2) 

MLE 

1.00765 

0.99696 

1.00145 

1.00132 

■  ■^■! 

( .559E-2) 

( .501E-2) 

( .488E-2) 

( .370E-2) 

CMLE 

1.00249 

1.00201 

1.00158 

0.99891 

( .556E-2) 

( .457E-2) 

( .400E-2) 

( .339E-2) 

CMLEBM 

0.96269 

1.00709 

0.99290 

0.99852 

(  .573E-2) 

(.458E-2) 

(  .428E-2) 

( .343E-2) 

FBLP 

1.00987 

0.99189 

1.00179 

1.00012 

( .858E-2) 

( .757E-2) 

( .689E-2) 

( .521E-2) 

EVLP 

1.00266 

1.00196 

0.99920 

1.00236 

( .846E-2) 

( .762E-2) 

(  .687E-2) 

( .538E-2) 
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a 

=  .1,  A  =  1 

.0 

Sample 

Methods 

Size\  20 

30 

40 

50 

CRLB 

1.00000 

1.00000 

1.00000 

1.00000 

( .500E-1) 

( .408E-1) 

(  .354E-1) 

( .316E-1) 

MLE 

0.98701 

0.97236 

1.00561 

1.01143 

>  ,0'^ 

( .588E-1) 

(  .510E-1) 

(  .491E-1) 

( .390E-1) 

CMLE 

1.01661 

1.00728 

1.01239 

1.01226 

( .561E-1) 

(  .453E-1) 

(  .409E-1) 

( .344E-1) 

CMLEBM 

1.06559 

1.02134 

1.01187 

1.00333 

( .568E-1) 

( .455E-1) 

( .411E-1) 

( .341E-1) 

FBLP 

1.02129 

1.03458 

1.00977 

1.01949 

•Hi 

( .784E-1) 

( .771E-1) 

( .680E-1) 

( .515E-1) 

EVLP 

0.98591 

0.97291 

0.99350 

1.00196 

W: 

( .800E-1) 

a 

( .751E-1) 

=  .5,  A  =  1 

( .679E-1) 

.0 

( .511E-1) 

Sample 

Methods 

Size\  20 

30 

40 

50 

CRLB 

1.00000 

1.00000 

1.00000 

1.00000 

(.250E-0) 

( .204E-0) 

( .177E-0) 

( .158E-1) 

MLE 

1.08614 

0.94289 

0.95671 

1.07197 

( .318E-0) 

( .287E-0) 

( .255E-0) 

( .243E-0) 

-*•  -r'l 

CMLE 

1.06013 

1.01013 

1.01601 

1.01038 

•;•>'■  5 

( .292E-0) 

( .218E-0) 

( .201E-0) 

( .169E-0) 

CMLEBM 

1.01685 

0.97094 

0.99316 

0.99043 

( .309E-0) 

( .234E-0) 

( .214E-0) 

( .177E-0) 

FBLP 

0.89376 

1.21643 

0.90185 

1.10967 

( .558E-0) 

( .534E-0) 

( .499E-0) 

( .320E-0) 

EVLP 

0.88750 

0.89921 

0.89958 

1.00196 

..-X;: 

( .562E-0) 

( .501E-0) 

( .468E-0) 

( .355E-0) 
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Table  6a 


Model  2 

Moving  Average  Error 

a  =  .1,  (ji  —  2 .0 


Sample 

Methods 

Size\  20 

30 

40 

50 

CRLB 

2.00000 
( .353E-2) 

2.00000 
( .193E-2) 

2.00000 
( .125E-2) 

2.00000 
(  .895E-3) 

MLE 

1.99998 
( .399E-2) 

2.00008 
( .279E-2) 

1.99999 
( .177E-2) 

2.00006 
( .171E-2) 

CMLE 

2.00015 
( .372E-2) 

2.00009 
( .201E-2) 

2.00001 
( .132E-2) 

1.99999 
(  .901E-3) 

CMLEBM 

2.00018 
( .381E-2) 

2.00000 
( .203E-2) 

2.00006 
( .140E-2) 

2.00001 
( .911E-2) 

FBLP 

1.99619 
( .757E-2) 

1.99658 
( .652E-2) 

1.99730 
( .427E-2) 

1.99641 
( .299E-2) 

EVLP 

2.00610 
( .779E-2) 

2.00186 
( .645E-2) 

1.99851 

(.>441E-2) 

2.00181 
( .315E-2) 

a 

Table  6b 
=  .1,  A  =  1. 

0 

Sample 

Methods 

CRLB 

Size\  20 

1.00000 
(  .460E-1) 

30 

1.00000 
( .372E-1) 

40 

1.00000 
( .333E-1) 

50 

1.00000 
( .289E-1) 

MLE 

0.96120 
(  .571E-1) 

1.00961 
( .511E-1) 

0.99280 
( .451E-1) 

0.99549 
( .377E-1) 

CMLE 

0.99256 
( .501E-1) 

1.00113 
(  .409E-1) 

1.00147 
( .367E-1) 

0.99990 
( .291E-1) 

CMLEBM 

1.00825 
( .517E-1) 

1.01048 
( .411E-1) 

1.00290 
( .375E-1) 

1.00034 
( .301E-1) 

FBLP 

1.03703 
( .827E-1) 

1.04686 
( .811E-1) 

1.05334 
( .723E-1) 

0.93356 
(  .700E-1) 

EVLP 

0.94707 
( .901E-1) 

0.93546 
( .859E-1) 

1 . 08556 
( .801E-1) 

0.91957 
(  .733E-1) 
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6. 


CONCLUSIONS : 


In  this  paper  we  try  to  estimate  the  parameters  of 
a  complex  sinusoidal  signal.  Since  any  algorithm 
analyzing  a  complex  signal  can  be  used  easily  for 
analyzing  the  corresponding  real  signal  (Kumaresan  and 
Tufts;  1982)  therefore  our  algorithm  can  also  be  used 
easily  to  analyze  the  corresponding  real  signal.  It 

is  well  known  that  the  estimators  depend  on  the  root  of 
a  certain  polynomial  equation.  It  is  observed  that  the 
coefficients  of  the  polynomial  exhibit  a  certain 
symmetry.  In  this  paper  we  modify  the  maximum 
likelihood  method  of  Kundu  (1993a)  by  using  this 
symmetry  constraints.  Bresler  and  Macovski  (1986)  and 
Kumaresan,  Scharf  and  Shaw(1986)  also  obtained  the 
maximum  likelihood  estimators  using  this  symmetry. 
Their  methods  are  very  much  similar  to  each  other  (see 
Bresler  and  Macovski;  1986).  Kumaresan,  Scharf  and 
Shaw  (1986)  suggested  to  use  some  standard  constrained 
optimization  packages  whereas  Bresler  and  Macovski 
(1986)  obtained  the  explicit  normal  equations  to  solve 
this  constrained  problem.  Since  both  the  methods  are 
almost  same,  we  report  the  results  of  Bresler  and 
Macovski  (1986)  only. 

The  numerical  results  confirm  the  satisfactory 
performance  of  the  CMLE  algorithm.  The  mean  estimators 
over  one  thousand  replications  show  that  the  sinusoidal 
frequencies  can  be  estimated  unbiasedly  by  all  the 
methods.  But  as  far  as  MSE  is  concerned,  the  CMLE  has 
the  lowest  MSE  in  almost  all  the  cases.  The  MSE  of 
the  CMLE  reaches  the  CRLB  in  some  situations.  As  N 

increases  or  a  decreases ,  the  MSE  of  the  estimators 
decreases.  When  N  is  large  and  a  is  small,  MLE,  CMLE 
and  CMLEBM  behave’ quite  similarly  but  for  small  N  and 
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large  value  a  CMLE  behaves  marginally  better  than 
CMLEBM.  It  is  observed  that  using  the  same  number  of 
iterations,  the  CMLE  has  a  clear  advantage  over  CMLEBM 
in  terms  of  MSE. 

Now  to  compare  the  computational  merit  of 
different  iterative  procedures,  observe  that  the  main 
computation  involves  in  each  iteration  in  inverting  a 
N-M  X  N-M  matrix  X*X,  where  X  is  of  the  form  (4.5). 
Since  X*X  is  a  banded  Toeplitz  matrix  with  bandwidth 
2M+1  the  inverse  computation  can  be  done  very 
efficiently  (see  Kumar; 1985  or  Kumaresan,  Scharf  and 
Shaw  (1986)).  In  each  iteration  to  calculate  CMLEBM  it 
is  reguired  to  solve  a  set  of  M+1  linear  equations  of 
M+1  unknowns  and  to  compute  CMLE  it  is  required  to 
obtain  the  eigenvalues  and  eigenvectors  of  a  real 
symmetric  matrix  of  order  M+1  x  M+1 .  Observe  that  in 
MLE  computation  we  need  to  calculate  the  eigenvalues 
and  eigenvectors  of  a  Hermitian ---matrix  of  order 
M+lxM+1.  So  CMLE  involves  less  computation  than  MLE 
and  CMLEBM  involves  less  amount  of  computation  than 
CMLE  in  each  iteration.  But  usually  M  is  very  small 
compared  to  N  in  practice,  so  computationally  it  does 
not  make  too  much  of  a  difference  in  each  iteration 
calculation. 

Another  important  problem  needs  discussion  is  the 
choice  of  initial  value  for  the  iterative  procedure.  It 
is  well  known  that  in  this  problem  unless  the  starting 
values  are  reasonably  good  it  is  very  likely  that  any 
iterative  procedure  converges  to  a  local  minimum  rather 
than  a  global  minimum.  The  work  of  Rice  and  Rosenblatt 
(1988)  showed  that  there  are  local  minimum  of  the  sums 
of  squares  with  respect  to  the  co  which  are  0(1/N) 

apart.  Since  the  EVLP  method  gives  a  convergence  rate 
of  only  0(l/yN),  one  would  not  expect  it  to  be  good 
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n 


enough  as  a  starting  value  to  find  the  global  minimum 
to  the  sums  of  sguares .  Although  we  observed  in  our 
simulation  study  that  the  EVLP  estimator  is  a 
reasonable  starting  point.  Bresler  and  Macovski  (1986) 
and  Kumaresan,  Scharf  and  Shaw  (1986)  also  used  the 
same  starting  point  and  obtained  global  convergence  in 
their  simulations.  We  observe  that  FBLP  also  can  be 
used  as  a  reasonable  starting  value.  We  perform  the 
experiments  with  both  the  starting  values  and  it  is 
observed  from  any  one  of  them  all  the  algorithm  namely 
MLE,  CMLE  or  CMLEBM  converge  to  the  same  global 
minimum.  Only  difference  it  makes  in  the  number  of 
iteration  and  that  is  also  very  marginal.  It  seems  if 
a  is  not  too  large,  then  EVLP  or  FBLP  estimator  can  be 
used  as  a  initial  guess  to  obtain  the  global  optimum. 
For  large  a,  it  is  not  very  easy  to  choose  a  good 
starting  value  without  having  any  prior  information 
about  the  parameters.  It  seems  more  work  is  needed  in 
this  direction. 

Observe  that  the  symmetry  constraint  of  the 
polynomial  coefficients  obtained  in  Section  3,  is 
necessary  but  not  sufficient  for  the  polynomial  roots 
to  be  on  the  unit  circle  (see  Kay  and  Marple;  1981). 
The  constrained  algorithm  may  return  real  roots  that 
are  reciprocal  of  one  another .  But  in  a  neighborhood 
of  the  true  value,  the  symmetry  constraint  is 
sufficient  also.  We  did  not  encounter  this  problem  in 
our  simulation  study  in  any  model. 

Recently  Kahn  et  al  (1992)  has  studied  extensively 
the  effect  of  different  scaling  factor  such  as  |g|  =  1 

or  q  =  1  etc.  on  the  consistency  of  Prony^s  method.  It 
o 

is  observed  that  different  scaling  factors  play  an 
important  role  in  the  asymptotic  inference.  It  can  be 
exploited  here  also.  It  seems  extensive  study  of 
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different  scaling  factors  is  needed  for  this  problem. 

Observe  that  the  MLE,  CMLE  and  CMLEBM  are 
asymptotically  eguivalent.  It  is  known  (Mitra  and 
Kundu;  1993,  Rao  and  Zhao;  1993)  that  the  asymptotic 
variance  of  the  maximum  likelihood  estimates  are  equal 
to  the  CRLB.  It  is  interesting  to  note  that  the 
variances  for  the  least  squares  estimates  are  rather 
greater  than  the  corresponding  CRLB  in  some  cases. 
Since  it  is  a  highly  non-linear  problem,  we  believe  it 
is  due  to  the  poor  asymptotic  approximation. 

Rice  and  Rosenblatt  (1988)  observed  that  unless 
the  frequencies  are  estimated  with  the  precision 
0(n“^),  the  corresponding  estimators  of  the  linear 
parameters,  say  by  (4.3),  may  not  be  consistent.  So 
although  MLE,  CMLE  or  CMLEBM  are  computationally  more 
expensive  than  FBLP  or  EVLP,  it  may  be  advisable  to 
use  them  to  obtain  consistent  estimates  of  both  linear 
as  well  as  non  linear  parameters.  Among  the  MLE,  CMLE 
and  CMLEBM  we  recommend  to  use  CMLE  one  since  it  has 
the  lowest  mean  squared  errors. 

APPENDIX 

In  this  Appendix  we  try  to  give  some  theoretical 
support  for  the  convergence  of  the  iterative  procedure 
described  in  (4.20),  particularly  when  the  errors  are 
independent  and  identically  distributed.  To  prove  the 
asymptotic  stability  of  the  iterative  process  (4.20), 
we  need  to  make  the  assumption  that  the  true  parameter 
value  (a,  <»>)  =  (o  ,•••/«  /W. /♦•  ww  )  is  an  interior 
point  of  the  parameter  space.  We  take  M  is  odd  for 
brevity.  The  proof  for  M  even  is  exactly  the  same. 
The  iterative  procedure  (4.20)  can  be  expressed  as 

F(c\  d*')  =  (c^'^S  d'^'^)  (A.l) 
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where  and  d*'  are  the 


first 


elements  of  the  vector  x  ,  i .  e .  x 
k 


and 

k 


the  last 


/  K 

=  (c  , 


k 

C  = 


(Cp'  *  *  • 


and 


d’‘=  (d’',...,d’' 


(M+l)/2 

where 

). 


Here  the  function  F:  K 


,M+1 


K 


,M+1 


is  defined  implicitly 


through  (4.20).  Therefore  the  convergence  matrix  F  of 
the  iterative  process  can  be  written  as 


(A. 2) 


(c,d)=(c,d) 

Here  x  =  (c,  d)  is  the  solution  of  (4.19)  and  all  the 
four  sub  matrices  are  of  the  order  (M+l)/2  x  (M+l)/2. 
The  (i/j)  th  element  of  dc'*^^^/dc'^  is  obtained  by  taking 

the  partial  derivative  of  with  respect  to  o’?  and 

the  other  elements  are  also  similarly  defined.  Finally 

the  matrix  F(c,  d)  is  obtained  by  evaluating  all  the 

elements  at  the  point  x  =  (c,  d) .  The  sufficient 

condition  that  x  is  a  point  of  attraction  of  the 
Iterative  process  (4.20)  is  that  the  spectral  density 

of  F(c,  d)  is  less  than  one  (Ortega  and  Rheinboldt, 
1970).  Writing  (4.20)  explicitly  using  the  expressions 


of  D  from  (4. 

11)  to  (4.14) 

,  we  obtain 

E  A.  .  c.  t 

‘  ^  ^-5  J 

J  =o 

<M-l>/2 

^  ^  jk+1 

E  B.  .  d. 

J  =o 

.  <k+l>  k+l 

=  X  c. 

(A. 3a) 

_  k+l  , 

E  B.  .  C.  + 

'  ^  J 

J  =o 

<M-l>/2 

E  c.  .  d. 

.  „  I-J  4 

4  =0 

+ 

11 

(A. 3b) 

1  ”  0^  1l  f  ••  9  f 

(M-l)/2  and 

=  1. 

Here  all  the  matrix  A/ B  and  C  are  evaluated  at  x  = 


F(c,  d)  - 


dc 

dc 


k+i 


dc 

dd^ 


k+l 


ddL 

^  J< 
dc 


k+l 


dd 

dd** 


k+l 
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(c’',  d^^).  The  elements  of  the  matrix  PCc,  d)  can  be 
obtained  by  differentiating  both  sides  of  (A. 3a)  and 
(A. 3b)  with  respect  to  c^  for  1  =  0,1,  .  .  .  ,  (M-l)/2  and 

A 

evaluating  at  the  point  x  yields 


A(c, 

Apj, 

d)  ^  +  B(c,  d) 

Arf* 

^  +  R,(c,  d)  = 

dc 

i 

dci 

== 

A  Afji 

dX  ^  dc 

dc  dc 

L 

0  > 

A,jt 

d)  ^  +  C(c,  d) 

+  S  (c,  d)  = 

dc^ 

dc^ 

= 

A  Afjt 

dX  ^  dd 

—  d  +  X  — 

dd,  dc 

1*  I* 

where 

R^(c,  d)  is  an  (M+l)/2  vector  whose  jth 

for  j 

=  0,  1, . . . , (M-l)/2  is  given  by 

(A. 4a) 


(A. 4a) 

element 


(c,  d)  = 

-  Y*U  (X*X)”^(UyX  +  X*U  )  (X*X)“‘X*Y 

J  I  L 

-  Y*U  fX*X)“^(u’'x  +  X*U.)  (X*X)"^X*Y 

u  J  J 

+  Y*x(x*x)"^(u[x  +  X*U^)(X*X)"^(U^X 
+  X*U.)(X*X)‘^X*Y  +  Y*X(X*X)“‘(u'^X 

J  ^ 

+  x*u.)(x*x)"‘(u7x  +  x*u  )  (X*X)"‘x*Y 

J  II 

-  Y*X(X*X)"^(u'^X  +  X*UJ  (X*X)'‘(U^Y 

-  Y*X(X*X)"^(UyX  +  X*U  )  (X*X)“^(U'^Y  (A. 5) 

U  U  J 

and  the  (M+l)/2  vector  S  (c,  d)  is  also  similarly 
defined.  Here  the  ith.  element  of  dc  /dc^  for  i  = 
0, . . . , (M-l)/2  is  obtained  by  taking  the  partial 
derivative  of  with  respect  to  c*'  and  evaluating  it 

at  (c,  d).  dd^'/dc  and  dx’^/dc  are  also  similarly 

t  ^ 
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defined. 

It  is  easy  to  see  that  because  all  the  quantities 
inside  the  parenthesis  in 


Lim  - 

f 

A 

+  \ 

^'"'1=0 

n  - >00 

L  dc^ 

dc  J 

L 

for  1  = 

0,1,  . 

•  •  / 

(M-l)/2 

(A. 6a) 

Lim  \ 

f  dX 
—  d 

A 

+  X 

o 

!1 

n  — >00 

dc 

i 

dc  ^ 
i 

for  1  = 

0,1,  • 

•  •  / 

(M-l)/2 

(A. 6b) 

both  (A. 6a)  and  (A. 6b)  are  bounded  over  a  compact  set. 
Therefore  from  (A. 4)  we  obtain 


A  A 

A ( c / d) 

A  A  A 

B*(c,d) 


A  A 

B(c,d) 

A  A 

C(c,d) 


dc 

dc 

dd 


A  A 

R(c,d) 

A  A 

_  S(c,d) _ 

>• 

=  0 

The  main  idea  of  the  proof  of  the  asymptotic 
stability  of  the  iterative  process  is  as  follows; 
(l/n)R^(c,d)  and  (l/n)Sj^(c,d)  converges  to  zero  vectors 
almost  surely  for  1  =  0,1, . . . , (M-l)/2  and  (l/n)D(c,d), 
where 


D(c,d)  = 


A( c , d) 

^  A  A 

B*(c,d) 


A  A 


B(c,d) 

A  A 

C(c,d) 


(A. 8) 


converges  to  a  positive  semi  definite  matrix  with  null 
space  spanned  by  (c*^,d°),  the  true  value  of  (c,d), 
corresponding  to  the  true  parameter  value  co.  Since 
cc’’  +  dd’’  =  1  (^*5) 

therefore 
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C  +  d  =  0  for  1  =  0,1,  ...,  (M-l)/2.  (A. 10) 

dc  dc, 

L  I 

Since  the  MLE  of  co  converges  almost  surely  to  the  true 
parameter  value  w,  therefore  (c,d)  converges  almost 
surely  to  (c°,d°).  From  (A. 10)  it  follows  that  the 


vector  [dc/dc  ,dd/dc  converges  to  a  vector  which  is 

L  1* 

orthogonal  to  (c°,d°).  Since  Lim  [dc/dc^ ,dd/dCj^ 

n  - >00  ,,  ,, 

belongs  to  the  null  space  of  lim(l/n)D(c,d)  and  the 

n-^00 

lateir  convBirges  to  a  positivG  sGiui  dGfinitG  matrix  with 

null  spacG  spannGd  by  (c^ ,d°)  ,thGrGforG[dc/dc^,dd/dc^] 
convGrgGs  to  zgito  VGctor  almost  surGly .  ThG  samG 


roasoning  holds  for  [dc/dd^,dd/dd^]*^  for  1  —  0,  1,..., 

(M-l)/2.  ThGrGforG  F,  thG  convorgonce  matrix, 
convGrgGS  to  ZGro  almost  suroly .  This  impliGS  thG 
asymptotic  stability  of  thG  proposGd  algorithm- 

To  complGtG  thG  proof  wg  nGGd  thG  following 
rGSults 


lim  ^  D(c,d)  =  lim  ^  E(D(c°,d°))  =  D 

— >00  ^  ^  ^  n — >-00 

lim  ^  R(c,d)  =  lim  i  E(R^(c  ,d  ) )  =0 
— >*00  ^  ^  ^  n-^00 

lim  ^  S^(c,d)  =  lim  ^  E(S^^(c  ,d  ))  =0 


n — >00 
for  1=0,  1, 


n — >00 
(M-l)/2 


Here  E  stands  for  expectation  and  is  a 


(A. 11) 
positive 


semi  definite  matrix  with  null  space  spanned  by  (c  ,d  ). 
The  proof  of  (A. 11)  is  very  much  similar  to  that  of 
Theorem  A.l  of  Kundu  (1994)  so  it  is  omitted. 
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